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ABSTRACT 

High-redshift clusters of galaxies show an over-abundance of spirals by a factor of 2 — 3, 
and the corresponding under-abundance of SO galaxies, relative to the nearby clusters. This 
morphological evolution can be explained by tidal interactions with neighboring galaxies and 
with the hierarchically growing cluster halo. The efficiency of tidal interactions depends on the 
size and structure of the cluster, as well as on the epoch of its formation. I simulate the formation 
and evolution of Virgo- type clusters in three cosmologies: a critical density model ^Iq — 1, an 
open model Hq — 0.4, and a flat model flo — 0.4 with a cosmological constant. The orbits of 
identified halos are traced with a high temporal resolution 10^ yr). Halos with low relative 
velocities merge only shortly after entering the cluster; after virialization mergers are suppressed. 

The dynamical evolution of galaxies is determined by the tidal field along their trajectories. 
The maxima of the tidal force do not always correspond to closest approach to the cluster center. 
They are produced to a large extent by the local density structures, such as the massive galaxies 
and the unvirialized remnants of infalling groups of galaxies. Collisions of galaxies are intensified 
by the substructure, with about 10 encounters within 10 kpc per galaxy in the Hubble time. 
These very close encounters add an important amount (10-50%) of the total heating rate. The 
integrated effect of tidal interactions is insufficient to transform a spiral galaxy into an elliptical, 
but can produce an SO galaxy. Overall, tidal heating is stronger in the low Qq clusters. 

Subject headings: galaxies: clusters: general — galaxies: evolution — galaxies: interactions 



1. Introduction 

This work is motivated by the wealth of recent 
observational discoveries of the strong morpholog- 
ical, chemical, and dynamical evolution of galaxies 
in clusters. Transformations of luminous galaxies, 
the existence of separate massive dark objects, and 
the evolution of cluster substructure intensify the 
need for theoretical understanding. Several differ- 
ent phenomena have been proposed as the "likely" 
explanation (see §2). While all of them take place 
in one form or another, it has eluded our under- 
standing which process is primarily responsible for 
the evolution. 

I argue that tidal interactions with neighbor- 
ing galaxies and with the dark cluster halo are the 
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favored scenario. Tides operate even before the 
cluster dynamically relaxes and affect both stars 
and gas in the infalling galaxies. Tidal effects add 
enough random motion to transform thin galactic 
disks into the kinematically hot thick-disk configu- 
rations characteristic of SO galaxies. Tidal disrup- 
tion of low density galaxies may provide a large 
fraction of the diffuse intracluster light. 

The relative importance of tides depends on 
the cosmological context. In an fio = 1 uni- 
verse, gravitating structures assemble late and the 
tidal effects are only now starting to play a role in 
shaping the galaxies, whereas in an open universe, 
ilo ~ 0-4, galaxies have formed early and experi- 
enced a substantial metamorphosis. The currently 
fashionable model, a flat universe with fl^ = 0.4 
and a cosmological constant, provides the inter- 
mediate case. In order to incorporate these effects 
into the galactic dynamics self-consistently, I per- 
form numerical simulations of the evolution of a 
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medium-size cluster of galaxies in the three cos- 
mologies. 

The choice of a numerical scheme is very im- 
portant here, as the required dynamical range of 
simulation is huge: from hundred parsecs to model 
the galactic dynamics to tens of Mpc to derive the 
cluster tidal field correctly. In order to achieve 
such a range with the existing codes, I introduce 
a new scheme of separating the cluster tidal field 
from the galactic model. First, I run a lower res- 
olution cosmological simulation that captures the 
details of cluster formation and its time-varying 
potential. At some early rcdshift, Zg = 5, I iden- 
tiiy the probable galactic halos and follow their 
trajectories until the present. Along each trajec- 
tory, I calculate and store the components of the 
tidal shear tensor. Then, I resimulate individual 
galaxies with a very high resolution applying the 
cluster tidal field as an external perturbation. The 
results of the stellar dynamical galaxy simulations 
are presented in a companion paper (Gnedin 2003, 
hereafter Paper II) . No single simulation of a clus- 
ter of galaxies (§2.5) has been able to achieve com- 
parable resolution on the galactic scale. 

I describe a careful procedure to extract the 
tidal field from the grid PM simulation. Because 
of the particle noise, it is necessary to apply a 
smoothing filter extending over a range of several 
grid cells. Also, large halos contribute significantly 
to the tidal field around their centers. Since we are 
interested only in the external tidal force on the 
galaxy, this self-contribution needs to be removed. 
Special care is taken to evaluate the contribution 
of the galactic particles and subtract it from the 
total tidal tensor. 

Also, the resulting tidal field depends on the 
correct calculation of the galactic trajectories. I 
use two independent methods to follow the tra- 
jectories: in one, the galactic centers of mass are 
assumed to be test particles following the overall 
cluster potential; in the other, the center of mass 
is calculated using the densest subset of particles 
constituting an actual halo. The sampling of halos 
at each time step provides an excellent time reso- 
lution of the tidal field, of order 10^ yr. Finally, 
the distribution of halos and their relative veloci- 
ties enable us to estimate the contribution of close 
encounters to the external tidal force. 

In §2, I review the present observational status 
of galaxies in clusters at various redshifts and the 



theoretical mechanisms connecting pieces of the 
puzzle. In §3, I describe the numerical method 
used to simulate the self-consistent evolution of 
the three clusters starting with the primordial fluc- 
tuations. In §4, I identify the galactic halos and 
calculate their trajectories. In §5, 1 study the pos- 
sible mergers and close encounters between galax- 
ies, including the effects of halo clustering. Finally 
in §6, I calculate the tidal field along the galaxy 
trajectories and use it in §7 to estimate the amount 
of tidal heating of galaxies. 

2. Scenarios of Galactic Evolution 

The idea that the environment plays a role in 
the evolution of galaxies has been recognized long 
ago. Spitzer & Baade (1951) suggested that col- 
lisions between galaxies might be responsible for 
transforming spirals into the gas-free SOs. Toomre 
(1977) advocated a merger of two spirals to pro- 
duce an elliptical. Richstone (1976) argued for the 
truncation of dark matter halos of galaxies by the 
tidal field of the cluster. Gunn & Gott (1972) 
proposed the stripping of gas by the ram pres- 
sure of the intracluster medium. Finally, Merritt 
(1983) and most recently Moore at al. (1998,1999) 
considered tidal interactions between the galaxies 
and the dense cluster core and between the galax- 
ies themselves. I review the observational picture 
and then discuss briefly the pros and cons of each 
of these scenarios. 

2.1. "Spirals vs Ellipticals" or "Spirals vs 

SOs"? 

The Hubble Space Telescope has been instru- 
mental in shaping our current understanding of 
distant clusters of galaxies. The unique resolving 
power of the HST brings us exquisite details of 
the distorted profiles of the otherwise normal disk 
galaxies. Strongly depleted in nearby clusters, spi- 
rals make up to 50% of the galaxies in clusters at 
redshift z ~ 0.5 (Oemler et al. 1997; Ellis et al. 
1997; Smail et al. 1997). These results are simi- 
lar to the Butcher-Oemler (1978) effect, an abun- 
dance of blue irregular galaxies at an intermedi- 
ate redshift. Spectroscopic evidence (Dressier et 
al. 1997) suggests bursts of recent, or even cur- 
rent, star formation. This picture is drastically 
different from the rich present-day clusters, popu- 
lated mostly by featureless elliptical or SO galaxies 
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which foUow the local morphology- density relation 
(Dressier 1980). It is well established now that a 
strong chemical and morphological evolution has 
taken place in clusters of galaxies. 

What are the options? Did the z = 0.5 spi- 
rals become the 2; = ellipticals? On the con- 
trary, the HST data show that the fraction of el- 
liptical galaxies increases with redshift, although 
not as strongly as that of spirals. Both of these 
effects occur at the expense of SOs. Their popu- 
lation is depleted by a factor of two or three at 
z ~ 0.5 (Dressier et al. 1997). Observations sug- 
gest that star formation has ceased in the ellipti- 
cals \sy z 0.6, but continued in the SOs at least 
until z ~ 0.3 (van Dokkum et al. 1998). The el- 
lipticals must have exhausted their supply of fresh 
gas at an earlier time, as indicated also by their 
very high central densities (Kormendy 1989). 

Even at low redshift the situation is far from 
being smooth (Conselice & Gallagher 1999). New 
ground-based imaging of the Coma and five other 
nearby clusters reveals a variety of distorted, 
tailed, "line", and dwarf galaxies. These arc indi- 
cators of the ongoing dynamical evolution. Many 
galaxies are disturbed and some show blue col- 
ors. Also, individual groups of dwarf galaxies 
have been found in Coma. These results were not 
expected for the smooth and regular clusters. 

A new observational puzzle has been presented 
by the discovery of massive dark objects in the 
cluster CL 0024-1-1654, through the technique 
of image reconstruction of gravitational Icnsing 
(Tyson, Kochanski & Dell'Antonio 1998; Kochan- 
ski, Dell'Antonio & Tyson 1998). The dark ob- 
jects arc not associated with any of the luminous 
galaxies in the cluster and may contain as much 
mass as 10^° Mq. One of the possible explanations 
is that these objects arc just "failed" galaxies, in 
which gaseous dissipation was insufhcient to form 
stars. Alternatively, they can be stripped dark 
matter halos of other galaxies. Analytical argu- 
ments and previous studies (§2.3) suggest that 
extended galactic halos may be effectively swept 
away by the tidal field of the cluster. 

2.2. Mergers of Galaxies 

Collisions of galaxies are assisted by the dense 

environment of the cluster. However, in a virial- 
ized cluster, galaxies move with the relative veloc- 



ities of the order 1000 km s^^ and even on a head- 
on trajectory will pass through each other without 
doing much harm to the stellar disks. The gas can 
be shocked and stripped away in such a collision 
(Spitzer & Baade 1951), which would leave two 
gas- free galaxies. But a simple estimate argues 
that close encounters are very infrequent (however 
see §5 and Appendix A). 

Consider a virialized cluster of galaxies with a 
one-dimensional velocity dispersion GcI = 1000 km 
s~^ and a virial radius Rd = 1 Mpc. Assuming 
that the cluster contains Ng = 1000 galaxies uni- 
formly distributed inside this radius, the number 
of close encounters a given galaxy of size i?g = 10 
kpc will experience in a Hubble time {tn = 10^° 
yr) is 

iVencl « TtRI V2acl tu « 1- (1) 

Here I have neglected gravitational focusing and 
assumed the relative velocity to be \/2aci- Thus a 

galaxy is expected to encounter one other galaxy 
over the course of its evolution. However, the 
probability of a merger of two galax;ies is much 
smaller: 

Pmer ~ A^enc,! « 10"^ (2) 

for the galactic velocity dispersion ag = 200 km s~' 
Here, three powers of Ogjad come from the phase 

space density and another power of Ug/ad comes 
from the limitation of the velocity that leads to a 
merger (Binney & Tremaine 1987). 

Also, ellipticals usually have much higher cen- 
tral densities than spiral galaxies. So if two spirals 
merge, it seems hard to close the gap in density 
without strong dissipation (Carlberg 1986). As 
Jim Gunn once said (Gunn 1987), "I do not think 
you can make rocks by merging clouds." Thus the 
merger scenario appears unlikely to dominate the 
dynamics of galaxies, although mergers could be 
important for the formation of the central giant cD 
galaxy (Ostrikcr & Hausman 1977; Merritt 1985; 
Lauer 1988) and the ultraluminous infrared galax- 
ies and AGNs (e.g., Canahzo & Stockton 2001). 

2.3. Truncation of Dark Matter Halos 

Richstone (1976), White (1976), and Merritt 
(1983), among others, have investigated the trun- 
cation of dark galactic halos by the tidal field 
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of the cluster. In the simplest case, where both 
the cluster and the galaxy halos have singular 
isothermal distributions, the tidal radius is ap- 
proximately 

Rf w — ijp, (3) 

where Rp is the distance of closest approach to the 
center of the cluster. 

Real clusters are not perfectly smooth. Any 
substructure left after the hierarchical formation 
of the cluster will amplify the truncation effect. 
It is expected to control the amount of dark mat- 
ter, and therefore, the mass of the galaxies. On 
the other hand, tidal truncation should not signif- 
icantly affect the dynamics of the galactic disks 
which arc much more compact. Morphological 
evolution of luminous galaxies must be due to 
something else. 

2.4. Ram-Pressure Stripping of Gas 

Gunn & Gott (1972) argued that the ram pres- 
sure of the intracluster medium (ICAI) might 
sweep away gas from the galactic disks. X-ray 
observations show the presence of a significant 
amount of hot gas at the virial temperature, con- 
tributing as much as 10% of the total cluster 
mass. While the stars move freely through this 
media, the gas is subject to the ram pressure 
-fram ~ PicmV , where picm is the density of the 
ICM and v ~ Ud is the velocity of the galaxy. In 
order to strip the gas, the ram pressure should 
exceed the restoring gravitational force per unit 
area of the disk: 

PicxnV^ > 27rGS,Sg, (4) 

where is the surface density of the stars, and 
Eg is the surface density of the gas. For exam- 
ple, if the gas with an average density 1 cm~^ is 
distributed over a scale height of 100 pc and the 
stellar density is typical of the solar neighborhood. 
Eg w 75 Mq pc~^, the ICM needs to be as dense 
as riicm > 4 X 10~^ cm^'\ Such densities are ob- 
served in the cores of rich clusters, and numerical 
experiments (Kundic et al. 1993; Abadi, Moore & 
Bower 1999: Quilis, Moore & Bower 2000) show 
that a typical spiral there may lose most of its 
atomic gas. 

Spiral galaxies in clusters have long been known 
to be deficient of atomic hydrogen (Krumm & 



Salpeter 1976; Haynes, Giovanelli & Chincarini 
1984). These "anemic" galaxies (van den Bergh 
1991) have normal molecular hydrogen abun- 
dances but lack the atomic component in the 
outer parts of their disks. Since the H2 gas re- 
sides in dense molecular clouds, it is more able 
to withstand the external pressure; on the other 
hand, the H i gas is usually distributed much more 
smoothly at lower densities and is susceptible to 
the ram pressure stripping. But Valluri & Jog 
(1991) noticed an interesting feature in three rich 
Abell clusters and the Virgo cluster: the H l defi- 
ciency is stronger for larger galaxies. Since bigger 
spirals (excluding the LSBs) usually have a higher 
surface density, they should be able to protect 
their gaseous component better. The observed 
trend is the opposite. Having considered various 
possible mechanisms of gas removal, Valluri and 
Jog conclude that the H i deficiency can only be 
explained by tidal interactions. 

2.5. Tidal Interactions 

Tides are a critical factor of the cluster envi- 
ronment. Tidal effects depend on relative sizes of 
the system and the perturber. Consider stellar 
systems from smallest to largest scales: stars are 
10^ times smaller than star clusters and experience 
almost no tidal effects (except in close binaries). 
Star clusters are 10^ times smaller than galaxies 
and experience significant tidal effects, which how- 
ever operate on a long time scale (Gnedin & Os- 
triker 1997). In contrast, galaxies are just 30 to 
100 times smaller than clusters of galaxies. There 
tidal effects ought to be important! 

Semi-analytic models of galactic encounters 
(Merritt 1983) predict the segregation of galaxies 
by mass and the formation of a central cD galaxy. 
The giant cD grows by swallowing other massive 
galaxies in a process that widens the gap between 
the first and the second ranking galaxies (Ostriker 
& Hausman 1977). Mass segregation in clusters 
is also supported by numerical simulations (Frenk 
et al. 1996). Recently, Dubinski (1998) confirmed 
that the central galaxy, rising by mergers of its 
neighbors, follows a de Vaucouleurs profile and 
would be classified as a giant elliptical. 

Restricted numerical simiilations of Byrd & 
Valtonen (1990) and Valluri (1993) show that tidal 
compression by the cluster halo produces spiral 
arms and tidal tails in some of the representative 
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models of disk galaxies. 

More recently, Moore at al. (1996b,1998) simu- 
lated the dynamics of galaxies in a fixed, isother- 
mal cluster. They found that fast close encoun- 
ters with massive galaxies tend to destroy many 
dwarfs, and even large galaxies lose their thin 
disks. So the spiral galaxies are essentially trans- 
formed into the ellipticals or dwarf spheroidals. 
However, the potential of a real cluster grows in 
time, which leads to two competing effects: (1) 
galaxies enter the cluster at a later epoch and have 
less time to interact (depending on cosmological 
model), and (2) the time- varying cluster poten- 
tial gives rise to large-scale perturbations in the 
tidal field and causes tidal shocks comparable to 
those from massive galaxies. New simulations by 
Moore at al. (1999) include the hierarchical clus- 
ter formation and show weaker evolution of galax- 
ies. Whether the galaxy-galaxy or the cluster- 
galaxy interactions dominate has not yet been un- 
derstood. 

Similar tidal effects operate on smaller scale in 

the Local Group and may lead to the morphologi- 
cal evolution of the dwarf satellite galaxies (Mayer 
et al. 2001). 

3. Numerical Method 

I use a Particle-Mesh code (Cen 1992) to sim- 
ulate the formation and evolution of three clus- 
ters of galaxies in the critical (fio = 1), open 
(Qo = 0.4), and flat (f2o + f^A = 1) cosmological 
models described in Table 1. The code solves Pois- 
son's equation using FFT on a rectangular grid 
and evaluates the density field using the Cloud-in- 
Cell (CIC) method. The code runs fast in paral- 
lel, allowing a large simulation with 256'^ particles 
distributed on a 512'^ grid. Some of the PM re- 
sults have been tested with the Adaptive Mesh 
Refinement code (Bryan et al. 1995). Specifically, 
I checked galaxy trajectories and the calculation 
of the tidal force. 

The comoving size of the simulation box is 
32 Mpc and the grid cell size is A = 62.5 
kpc. This is a compromise between the desire 
to reach a reasonably high resolution in the clus- 
ter center and the need to have enough compu- 
tational volume to include late infall. The lat- 
ter is important because a very strong initial den- 
sity perturbation may lead to more mass being 



bound to the cluster than is available in the whole 
box, 9.1 X 10^^ Qoh~^ Mq, causing the cluster to 
"under-form" . 

Tests with different amplitudes of the power 
spectrum, as, show that the dark matter veloc- 
ity dispersion artificially saturates for ug » 1 in 
the low Oq models. Therefore, I adopt cjs = 1 for 
these models, in rough agreement with the relation 
derived from the comparison of largo scale TV-body 
simulations with the all-sky survey of X-ray clus- 
ters, as w 0.5 (Eke, Cole & Frenk 1996; Pen 
1998). Then the virial mass of the simulated clus- 
ters is about 4 x 10"'^^ Mq, or 7% of the total mass 
within the computational box. The actual amount 
of mass bound to the clusters may be a factor of 
2-3 higher, up to 20% of the total mass. In the 
critical density model, the normalization ag = 0.5 
is chosen such that the cluster velocity dispersion 
is the same as in the other two models. 

The simulations start at redshift Zi = 20. The 
time step is determined by the requirement that 
no particle crosses more than a half of the grid 
cell, 0.5A. Also, it limits the expansion parameter 
a = {1 + z)~^ to grow by no more than 0.05(1 -|- 
Zi)~^. I have checked that increasing the time step 
by a factor of two does not change noticeably the 
structure of the resulting cluster. 

The look-back time as a function of redshift 
(Peebles 1993, p. 314) can be calculated analyt- 
ically in each of the three models. For fio = 1, 
the standard result is Hq t{z) = | (1 + z)~^/'^. For 
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In 
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/ 



f = 



I + Q-qz 
1-J^o ■ 



1-fio 



(6) 

The Hubble constant is assumed to be Hq = 65 
km Mpc"-^ {h = 0.65), according to the time 
delay measurements in a gravitational lens system 
0957-F561 (Kundic et al. 1997) and other indepen- 
dent methods (Freedman et al. 1994). 

On 16 processors of SGI Origin 2000, the full 
simulation with a 512^ grid runs for about 36 
hours. A large fraction of the time is spent on 
the calculation of the tidal field due to and around 
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the identified halos. All simulations presented in 
this paper were done using the 64-processor SGI 
Origin 2000 supercomputer at the Princeton Uni- 
versity Observatory. 

3.1. Constrained Initial Conditions 

The cluster simulations use the constrained 

Gaussian random initial conditions, generated us- 
ing the method of Hoffman & Ribak (1991). The 
actual realization is taken from the cluster com- 
parison project of Frcnk ct al. (1999). In that 
work the initial density perturbations were cho- 
sen to produce a 3a peak with a Gaussian filter 
of radius 5h~^ Mpc. Retaining the same random 
phases from the original realization, I have renor- 
malized the perturbations according to the power 
spectrum of the chosen cosmological models. 

For each of the three models, I evaluate the 
power spectrum P{k)nev/, including the effects of 
baryons, using the fitting expressions of Bardeen 
et al. (1986) and Hu & Sugiyama (1996). On the 
cluster and sub-cluster scales, the power spectra 
follow approximately a power law P{k) cx fc~^'^. 
Then I transform the initial displacements d into 
the Fourier space and renormalize their ampli- 
tudes according to the new power spectrum: 

jnew _ jold f P{k)ncw \ ^ f„\ 

where the old power spectrum, P(fc)oid: is from 
Prenk et al. (1999). Both power spectra are nor- 
malized to the corresponding value of <y1{zi) at the 
initial redshift Zi. Finally, I transform the per- 
turbation field back to real space and obtain the 
initial displacements for the dark matter particles. 

3.2. Structure of the Clusters 

The parameters of the resulting clusters are 
given in Table 2. The virialized region is taken 
to be a sphere originating at the cluster center, 
with a mean overdensity 5 = 200. The average 
mass of the clusters is 4 x 10^^ Mq, and the one- 
dimensional velocity dispersion is 660 km s~^. 

Figure 1 shows the projected surface density 
of the three clusters at the present epoch. The 
fio = 1 cluster has not yet relaxed and still contin- 
ues to accrete matter. Note a large infalling group 
of galaxies to the upper right from the center on 
the enlarged panel. In contrast, the two low Oq 



clusters appear more spherical and regular, and 
also very similar to each other; the cosmological 
constant does not seem to affect the cluster prop- 
erties. 

Figure 2 shows that the density profiles of the 
three clusters are very similar, with a slight excur- 
sion of the Q.Q = 1 cluster at r « 2/i~^ Mpc due 
to the infalling group. This group is clearly visible 
on the velocity dispersion plot, Figure 3. 

The decrease of the velocity dispersion at small 
distances from the center is partly real and partly 
numerical. The high resolution dark matter simu- 
lations in Prenk et al. (1999) show a similar de- 
crease; therefore our models should exhibit the 
trend. Also, such a decrease exists in the popular 
NFW model (Navarro, Frenk & White 1997). On 
the other hand, present simulations lack resolution 
in the cluster core and cannot be trusted at small 
radii. To illustrate the resolution effects, I have 
repeated the simulations of Cluster III on smaller 
grids, 256^ and 128^ cells (Table 3), with the corre- 
spondingly fewer particles, Ng^id = 2^ A^p. Figure 
4 shows that the density profile is reproduced well 
to the resolution limit of each simulation, but the 
velocity dispersion starts to decline much earlier in 
the lower resolution runs. Nevertheless, the differ- 
ence between the 512'^ and the 256"^ runs is small 
and confirms the accuracy of the larger simulation. 

The simulated clusters are similar to the Virgo 
cluster in size and velocity dispersion. ROSAT 
results (Bohringer 1994) indicate that Virgo's X- 
ray mass is (1.5 — 5.5) x 10^** M©, while the opti- 
cal velocity dispersion of galaxies is between 581 
and 721 km s~^ (Huchra 1985). These values are 
quite common for the nearby clusters, as shown 
by the optical and X-ray survey of Girardi et al. 
(1998). Often cited in theoretical work, the Coma 
cluster is much larger, with the one-dimensional 
velocity dispersion about 1000 km s~^ (Merritt 
& Gebhardt 1994) and the virial mass of order 
2 X 10^''' Mq (Bohringer 1994). Our models, there- 
fore, represent modest but more common clusters 
of galaxies. Their velocity dispersion is character- 
istic of the Abell richness class 1 (see Postman, 
Huchra & Geller 1992). 

4. Galaxies in the Simulation 

The definition of a galaxy in the cluster simu- 
lation is not trivial. While at high redshift dark 
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matter halos are significantly denser than the sur- 
rounding media, at low redshift the density con- 
trast effectively disappears. This effect is known 
as overmerging of halos (e.g., Moore, Katz & Lake 
1996a; Klypin et al. 1999) due to the insufficient 
mass resolution. Therefore, it is difficult and, 
to some extent, uncertain to try to identify the 
galaxy halos from the particle distribution at all 
times. Instead, I have implemented the follow- 
ing simple procedure, which can be extended to a 
higher degree of sophistication. 

I identify the dark matter halos at a time early 
enough that the density contrast is large, but late 
enough that galaxies are known to exist. A fiducial 
epoch at which the central parts of large galaxies 
are assumed to have formed is taken to be a;^ = 5. 

In order to obtain a statistically significant sam- 
ple, I choose 100 most massive halos identified 
by the group finding algorithm. Then, I add 100 
massless particles to the PM simulation with the 
coordinates and velocities of the halo centers of 
mass. These new "galactic" particles do not affect 
either the density or the potential calculation but 
rather act as test particles. I follow the motion 
of these particles and output the tidal field along 
their trajectories at each time step. 

By not using real particles to calculate galaxy 
trajectories, we effectively avoid the overmerging 
problem: the test particles will not merge even if 
they come very close to each other. The drawback 
is that due to a finite time step the test particles 
may gradually depart from the bulk of the galactic 
particles. To guard against this effect, I have also 
implemented another scheme. 

The second scheme involves real particles con- 
stituting the halo. I calculate the halo density at 
each particle's position using the same CIC kernel 
that is used in the main PM calculation. Then I 
take the densest subset of particles, 1 /8 of the to- 
tal, to calculate the center of mass directly. These 
dense particles presumably lie deep in the poten- 
tial well of the halo and serve as robust indicators 
of the center. In case of the two galaxies merging, 
their centers defined this way would clearly indi- 
cate it. However, I do not attempt to resolve the 
internal structure of galaxies in these simulations 
and only trace their centers of mass. 



4.1. Group Finding Algorithm 

Dark matter halos arc identified using the 
coordinate- free algorithm HOP proposed by Eisen- 
stein & Hut (1998). It is similar in spirit to DEN- 
MAX (Gelb & Bcrtschinger 1994) but uses densi- 
ties at the particle positions instead of the regular 
grid. I have extended HOP by including a pro- 
cedure for the removal of particles which are not 
gravitationally bound to the group. This proce- 
dure, similar to the SKID algorithm (Weinberg, 
Hernquist & Katz 1997), is outlined below. 

The gravitational energy of the particles in a 
group is calculated by direct summation of the 
Plummer potentials with a smoothing length cor- 
responding to one grid cell of the cluster simu- 
lation. The kinetic energy includes the Hubble 
expansion term. The total energy per unit mass 
is 

-Ej = - [Vi - Vcom + Hg{ri - rcom)] (-(-J.. _ i.^.)2!^ g2)l/2' 

(8) 

where m.p, r^, v^ are the particle mass, physical 
coordinate, and velocity, respectively, and Hg is 
the Hubble constant at redshift Zg = 5. The index 
"com" refers to the center of mass of the group. 
The particles are unbound as follows. I find par- 
ticle m with the maximum energy and discard it, 
if Em > 0. I recalculate the center of mass using 
the remaining particles and update the potential 
energy. Then, I find a new highest energy par- 
ticle and repeat the procedure. In the end, all 
remaining particles are bound or else the group is 
dispersed. 

A straightforward implementation of the un- 
binding procedure works well, although somewhat 
slowly for groups containing more than 10^ parti- 
cles. In the longest case, the unbinding algorithm 
took about 40 minutes on a single SGI RIOOOO 
processor. 

The second part of the HOP algorithm, RE- 
GROUP, uses three different overdensity param- 
eters to define the outer boundary of the group 
(Souter), the saddle point for merging two groups 
together {5 saddle), and the minimum density of a 
viable group {Speak)- After some experimenting, I 
have chosen the following parameters to produce 
the reasonable-size halos: [Souter, 5 saddle, Speak] = 
[15,45,50] for the Q.o = l model, and [20,55,60] for 
the low Oo models. The different sets of parame- 
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tcrs arc necessary because of the different masses 
of particles (wp = 5.4 x lO^fio/i"^ Mq). 

Figures 5-7 illustrate the location of particles in 
the identified halos. Notice the filamentary struc- 
ture of bound particles indicating the directions 
along which the halos continue to assemble. Most 
galaxies in the central part of the simulation vol- 
ume (shown on the Figures) fall eventually into 
the forming cluster. Because of the difference in 
particle masses, the fio = 1 halos contain signifi- 
cantly fewer particles than the = 0.4 halos. 

4.2. Halo Mass Function 

The mass function of the identified halos can 
be compared with the prediction of the Press- 
Schechter theory (Press & Schechter 1974). I have 
calculated the expected number density of galax- 
ies at Zg using equations (A4-A9) of Mo, Mao & 
White (1998). 

Figure 8 shows the halo mass functions in the 
three clusters. Note the relative depletion of low 
mass halos in the simulations. This is partly due 
to the truncation of the input power spectrum at 
the Nyquist frequency, which effectively prevents 
the formation of small-scale objects. For compari- 
son, the lower right panel shows the mass function 
obtained in the lower resolution, 256^ simulation 
of Cluster III. As expected, the smaller halos are 
depleted more strongly than the larger ones. 

4.3. Galaxy Trajectories 

Figures 9-11 show the trajectories of two se- 
lected halos in each cluster. These galaxies are 
used for the high resolution simulations in Pa- 
per II, representing a large spiral and a dwarf 
spheroidal, respectively. Their trajectories are 
quite irregular. In the f^o = 1 case, the galaxies 
are just entering the cluster at the present epoch, 
although not on a free-fall path. The orbits are 
similar to those shown in Frenk et al. (1996). On 
the other hand, in the low fio cases the galax- 
ies have enough time to move around the cluster 
center several times, more so in the open model. 
These differences will have major effects on the 
dynamical evolution of the galaxies. 

In general, there is a good agreement between 
the trajectories calculated using two different 
methods, especially for the less massive galax- 
ies. For the very massive halos, the "test particle" 



motion occasionally tends to oscillate around the 
trajectory determined by the densest particle sub- 
set. I have therefore decided to use the latter as 
the more reliable trajectories. 

5. Mergers and Close Encounters 

The hierarchical formation of clusters and the 
survival of substructure lead to more frequent col- 
lisions of galaxies and lower relative velocities at 
early epochs. Preferentially radial infall implies 
that the final distribution of galaxies is very dense. 
Most of the galaxies lie close to the cluster center, 
while the rest form a scarcely populated halo. 

The number density of massive galaxies (Mmin ~ 
10^^ Mq) can be fitted by a power law with a cen- 
tral core: 



ng{r) 



T7.0 



(9) 



(r + ro)T' 

normalized to the total number of galaxies Ng 
within the virial radius of the cluster, i?vir- Not 
all initially identified halos end up in the clus- 
ters, so with Ng < 100 trying to estimate the 
core radius ro and the slope 7 by fitting a his- 
togram of galaxy positions yields highly uncer- 
tain results. Instead, I use a maximum likelihood 
(ML) method analogous to that in Gnedin (1997). 
The parameters ro and 7 are to be determined 
by searching for the maximum of the ML func- 
tion, C = e~'^i> riiii 47rrfn(ri). Requiring the two 
partial derivatives of log/!(ro,7) to vanish gives 
a system of two integral equations, which have 
been solved numerically using the software pack- 
age Maple. The resulting parameters are given in 
Table 4. 

The slopes of the density profile are steep, 
7 = 2.4 — 4.2. The massive galaxies are more con- 
centrated in the inner parts of the clusters than 
the dark matter which follows an approximately 
isothermal profile, pdm r^^. Also, the core 
radii in all three clusters are very small, below 
the force resolution scale. At such small distances 
from the center the distribution of galaxies can be 
affected by the giant cD galaxy which may swallow 
its neighbors. To check this, I have repeated the 
ML analysis excluding all galaxies within one res- 
olution element. A, from the center. This forces 
the core radii to be at least one A. The resulting 
values (Table 4) are very close to this limit, sug- 
gesting that the true core radii can indeed be quite 
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small. Thus the effect of the central giant on the 
overall galaxy distribution seems to be small. 

5.1. Number of Close Encounters 

The detailed trajectories of galaxies can be used 
to calculate the number of close encounters during 
the epoch of cluster formation. In order to detect 
all encounters, the trajectories are linearly inter- 
polated between the successive time steps. The 
comoving separations are converted into physical 
kiloparsecs using the appropriate value of the Hub- 
ble constant at each time step. If a galaxy is found 
to have merged (see Appendix A), its subsequent 
trajectory is ignored. 

Table 5 shows the number of encounters, Ngnc, 
as a function of the impact parameter, i?enc- It is 
defined as a minimum distance between the cen- 
ters of mass of the galaxies. The last column gives 
the "total" number of encounters for i?enc = 40 
kpc (it is assumed that no mergers are possible 
for the separations larger than 40 kpc). Table 5 
also shows the average number of close encounters 
per galaxy, A^onci- Since each encounter involves 
two galaxies, this number is N^nc,! = "2 Ncnc/Ng. 

Figure 12 shows the number of collisions per 
galaxy within i?enc = 20 kpc. Note that in the low 
flo clusters, most collisions take place near the end 
of the simulation, when the galaxies have assem- 
bled inside the cluster. In Cluster I, the galaxies 
are closer to each other initially and have close 
encounters uniformly distributed in time. Several 
distinct peaks, notably in Clusters I and III, can 
be due to the new groups of galaxies entering the 
clusters. 

The number of close encounters does not scale 
very well with R^^^, as predicted by equation (1). 
Since most encounters take place in the already 
formed cluster, any mergers prior to that remove 
prospective galaxies and reduce the number of 
subsequent collisions. To test that equation (1) 
would have been valid in the absence of merg- 
ers, I count the number of encounters disallowing 
mergers. The new numbers, Nenc, scale as R^^^ in 
agreement with the prediction (see Table 5). 

However, the total number of encounters is 
much higher than expected. The effect is mainly 
due to the concentration of galaxies at the center 
of the cluster, implying a smaller effective volume. 



Equation (1) can be rewritten in a general form: 

-^enc — 

j nRlVenctHnl{r)dV, (10) 

where the integration extends over the virial vol- 
ume of the cluster, V. The mean value of the num- 
ber density squared, {n"^), is substantially higher 
than the square of the mean density, (rig)'^, be- 
cause of the galaxy clustering. Their ratio is the 
clumping factor 

Q — ' g' J 9^ ' fill 

Using the azimuthally-averaged density profile (9) 
and the auxiliary functions 

I find C = fof2/fi (see Table 4). Now, the number 
of close encounters is 

(13) 

This is much higher than the naive estimate (1) 
and agrees well with the results for Cluster HI. 

Because of the higher degree of clustering in the 
low Oq models, more than half of all encounters 
occur within the inner 0.4 Mpc from the center. 
In contrast, in the = 1 model collisions take 
place throughout the cluster, leading to a smaller 
value of the clumping factor. 

5.2. Relative Velocities 

Figure 13 shows the relative velocities of galax- 
ies at the time of collision. The distributions 
are skewed, peaking at low values and extending 
large tails into the high velocities. A characteristic 
value, the median, is around 800 km s~^. Again, 
there is a strong difference between Cluster I and 
the low Oo clusters. Table 6 demonstrates that 
while the peak velocity is higher in Cluster I, the 
median and mean values are substantially lower 
than those in Clusters II and III. The very high 
velocities are almost unattainable in the = 1 
model. 

Table 7 shows that the velocity dispersion of 
galaxies is smaller than that of dark matter and 
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varies strongly with the position in the cluster. 
Due to a small number of identified galaxies I 
have split them in three equally-populated radial 
bins. In Clusters II and III, the velocity gradi- 
ent is strong and inverted, with the lowest disper- 
sion within the inner cluster core. In Cluster I the 
effect is less conspicuous. Since the samples are 
small. Table 7 serves only to indicate the trend. In 
contrast, "test particle" galaxies tracing the dark 
matter potential show a normal gradient weakly 
rising towards the central bin (from 550 km s~^ 
to 750 km in Cluster II). Thus the inverted 
gradient is a property of only real galaxies. 

A review of the literature reveals that numeri- 
cal simulations often produce "galaxies" with the 
pairwise velocity dispersion lower than that of 
dark matter (Colin, Klypin & Kravtsov 2000; 
Summers, Davis & Evrard 1995; Carlberg 1994; 
but see also Okamoto & Habe 1999; Diaferio et 
al. 1999; Ghigna et al. 2000). In a specific study 
of the relative velocities at the time of galaxy col- 
lisions, Tormen, Diaferio & Syer (1998) find that 
the peak velocity is only Vrei ~ 2.5 ag « 500 km 
s~^ (their Fig. 7). There are several reason for 
this: (i) A statistical bias of hierarchical structure 
formation, whereas galaxies form around highly- 
clustered density peaks. Therefore, large galax- 
ies are always around the cluster core and do not 
gain velocity from gravitational infall. (ii) A phys- 
ical bias relating to gas cooling and star forma- 
tion (Cen & Ostriker 2000) makes younger galax- 
ies move more slowly than the older ones, (iii) Dy- 
namical friction reduces the angular momentum of 
large galaxies and draws them towards the center. 
In the present simulations, the final distances from 
the center are anti-correlated with galaxy masses 
(Spearman's correlation coefficient is Vs ~ —0.3 
for Clusters II and III, but only Vg « —0.1 for 
Cluster I). Thus there is a moderate amount of 
mass segregation in the clusters formed at earlier 
times. 

Finally, numerical effects may lead to low galac- 
tic velocities. The PM code loses force resolution 
in the two inner grid cells, within which the galax- 
ies do not "feel" each other's gravity. Also, if the 
halo particles departed from its center, the center- 
of-mass velocity would be significantly reduced. 
However, if they were scattered completely ran- 
domly in the cluster the average velocities would 
have fallen below 100 km s~^, which is clearly 



ruled out. 

At present, observations neither confirm nor re- 
ject the velocity bias. Carlberg, Yee & EUingson 
(1997), in a survey of 1150 galaxies from 14 clus- 
ters, find an almost constant projected velocity 
dispersion up to the virial radius. The line-of-sight 
projection may, however, wash out a decrease of 
the true velocity dispersion near the center. An- 
other large survey of 15 rich clusters (Dressier 
& Shectman 1988) shows a significant amount of 
substructure with the lower pairwise velocities of 
galaxies, as well as a modest decline (of the order 
20%) of the azimuthally averaged dispersion in the 
central regions. So the current observational sta- 
tus of the velocity bias is inconclusive. 

5.3. Merger History 

The low relative velocities of galaxies make 
mergers more likely. Without having enough reso- 
lution to study mergers in detail, I adopt a simple 
prescription based on the impact characteristics of 
the galaxies. It consists of the three merger crite- 
ria described in Appendix A. 

Figure 12 and Table 5 show the expected num- 
ber of mergers per galaxy, A^mer,i- It is calculated 
as N^e,-/Ng, since each merger removes only one 
galaxy. The calculated merger probability is high. 
At impact parameters i?cnc < 10 kpc it is already 
A^mer,! ~ 10%, and the total value reaches > 20%. 
Most of the mergers take place shortly after galax- 
ies enter the clusters (even earlier in Cluster I). 
Once the clusters virialize, galactic velocities be- 
come too high and prohibit mergers. Yet, the 
number of mergers is only a very small fraction, 
about 0.1%, of the number of close encounters (cf 
eq. [2]). 

Ghigna et al. (1998) find very few mergers 
within the virial radius of their high-resolution 
cluster. This can be reconciled with our results 
because Ghigna et al. (1998) count only the re- 
cent encounters in the virialized cluster (z < 0.5) 
while our numbers include all encounters starting 
from Zg = 5. Note that most mergers occur at rel- 
atively early times (Figure 12). In the three cases 
considered, the last merger occurred at 2.3, 1.0, 
and 3.3 Gyr before the present, respectively. The 
merger probability of 22-29% may reflect the in- 
stances of the formation of large elliptical galaxies 
at high redshifts. They are expected to take place 
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naturally in the centers of future clusters. 

6. Tidal Field Around Galaxies 

The tidal field of the cluster is calculated along 
the trajectories of the identified galaxies. The 
tidal force Ftid,a and the shear tensor Fajs are 

where $ is the external potential and r is the 
radius-vector in the galactic reference frame. By 
definition tidal force is zero at the center of mass 
of the galaxy. I evaluate the tidal tensor at each 
time step for each identified halo, using the tri- 
linear interpolation scheme with a smoothing filter 
described in Appendix B. 

Being coordinate-independent, the trace of the 
tidal tensor serves as a good measure of the tidal 
field. Figures 14-16 show the trace F^a aroimd 
the pairs of large and dwarf galaxies chosen in §4.3 
from each cluster simulation. The time variation 
of the tidal force is extremely irregular. There are 
many strong peaks, reaching the amplitude of 200 
Gyr~^ and higher. The largest peaks are usually 
narrow, with the characteristic duration of 10^ yr, 
while the smaller peaks last longer. The occur- 
rence of the two types of peaks is clearly differen- 
tiated in time: the larger happen early, when mas- 
sive galaxies are distinctly denser than the back- 
ground; the smaller happen at later stages, when 
the halos merge and the gravitational potential be- 
comes smoother. 

The amplitude of the tidal field is the largest in 

the Qq = 0.4 simulation, although strong pertur- 
bations continue for a longer time in the A model. 
In the Oo = 1 model, the overall amplitude is 
smaller but there is more temporal variation. The 
main distinction between the galaxies in the low 
and high Oq models is the time spent within the 
cluster. In the latter case, the galaxies arc just 
arriving to the center, while in the low models 
the galaxies have entered earlier and experienced 
more tidal interaction. 

The overall amplitude of the tidal force is de- 
termined by the epoch of cluster formation. For 
the same virial overdensity, the cluster is denser 
on average when the mean density of the Uni- 
verse is high. In the low flo models clusters form 
at a higher redshift and therefore exert stronger 



tidal forces than in the Oo = 1 model. Since the 
subsequent growth of the clusters is balanced by 
the drop of the mean density, the tidal force is 
expected to be of the same order of magnitude 
throughout the evolution. 

Note that the peaks of the tidal force are not 
always correlated with closest approach to the 
cluster center. The galaxies are being tidally 
"shocked" at 1 Mpc from the center as often as 
at 400 kpc. Even though the center of the clus- 
ter has been moving slowly in the simulation box, 
it is still likely that close encounters with massive 
galaxies or groups of galaxies happen relatively far 
from the center. As the cluster grows, such tidal 
encounters occur within a fairly large (comoving) 
volume. 

The amplitude of the tidal force tells us about 
the effective scale of perturbation. For a smooth 
cluster, the tidal tensor is of the order 

GMjr) a% 

^ 72-' (15) 

where M(r) is the encloscid mass within radius r. 
In the simulated clusters the three-dimensional ve- 
locity dispersion is ad = a/S x 660 km s~^. If a 
characteristic scale, the cluster core radius, is 200 
kpc then the expected tidal perturbation is 

^ 32 ( ,V ( ' Y Gyr-\ 

VllOOkms-V V200 kpcy' ^ 

(16) 

The measured amplitude of the tidal force (cf. Fig- 
ures 14-16) is about 150 Gyr~^. Turning the argu- 
ment around, such tidal force corresponds to the 
effective scale of perturbation 

r, = 92 f ( , \ kpc. 

^ VllOO km s"V V150 Gyr-V 

(17) 

Smaller than the core radius but larger than indi- 
vidual galaxies, this scale is appropriate for small 
groups of galaxies. As discussed in §2.1, groups 
maintain their identity in many clusters. Even 
such a relaxed cluster as Coma shows a detectable 
amount of substructure (CoUess & Dunn 1996; 
Conselice & Gallagher 1999). Since large peaks of 
the tidal force dominate the total amount of heat- 
ing (see §7), I conclude that the tidal effects are 
produced primarily by the interaction of galaxies 
with the remaining substructure in clusters. 
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7. Tidal Heating 

The main effect of external tidal forces on a 
galaxy is a positive energy change of stars and 
dark matter, or tidal heating. The amount of tidal 
heating can be estimated analytically, at least for 
the spherically-symmetric component. For a sin- 
gle tidal shock, the energy input is the ensemble 
average of the velocity change squared, |((At;)^), 
where Av is a time integral of the tidal force times 
the reduction factor that takes into account the 
conservation of adiabatic invariants of stellar or- 
bits (Gnedin, Hernquist & Ostriker 1999). In the 
simulations, the tidal force has many peaks of vari- 
able duration. In addition, different components 
of the tidal tensor are not perfectly correlated so 
that they experience maxima and minima at dif- 
ferent times. Clearly, this complex behavior can- 
not be fully described as simply as a single tidal 
perturbation. Instead, I design a cumulative pa- 
rameter, based on the semi-analytic theory of tidal 
shocks, which provides a useful scaling relation for 
the amount of tidal heating. 

Assume that each peak of the tidal force can 
be considered a separate "tidal shock" . The total 
amount of heating is then the sum of the contribu- 
tions from all peaks. Moreover, since the velocity 
changes add in quadrature, all components of the 
tidal tensor can be treated separately. Then it fol- 
lows from the semi-analytic theory (cf. Gnedin & 
Ostriker 1999) that the total energy per unit mass 
of the system changes by 

(A£;) = l/tid(r^), (18) 

where r is a characteristic size of the galaxy and 
/tid is the tidal heating parameter 

n a,l3 / n \ dyn/ 

(19) 

where the sum extends over all peaks n and all 
components of the tidal tensor, a,/3 = {x,y,z}. 
Here t„ is the effective duration of peak n for each 
value of a and /3, and tdyn is the half-mass dynam- 
ical time. For the galaxy model from Appendix A, 
the dynamical time scales with the halo mass as 
tdyn C)C Mg^*. The inner integration is limited to 
the time intervals between the successive minima 
of the tidal force. In carrying out the integration. 



I have linearly interpolated F^fj between the out- 
put points to match the time step of the galactic 
simulations presented in Paper II. This provides 
the most direct way to compare the analytical esti- 
mate /tid with the results of self-consistent TV-body 
simulations. The agreement is generally good. 

Figure 17 shows the distribution of parame- 
ters /tid for the three clusters. There are more 
galaxies with the low value of /tid than with the 
high - tidal effects are different for different galax- 
ies. The histograms display a distinct trend: the 
low VLq clusters have stronger tidal heating than 
the fio = 1 cluster. A least-squares fit to the 
histograms, using a convenient functional form 
-/V(/tid) = No exp(-/tid//o), yields the following 
characteristic parameters: 

/o = 8 X 102 Gyr-2 (Cluster I), 

/o = 5 X 10^ Gyr-2 (Cluster II), (20) 

/o = 2 X 103 Gyr-2 (Cluster III). 

The timescalc of perturbation, t„ , is very impor- 
tant in determining the total amount of heating. 
Without the adiabatic correction factor, the values 
of /q would have been at least an order of magni- 
tude higher. Overall, the open Q.q = 0.4 model has 
the highest heating rates and the Oq = 1 model the 
lowest. Thus galaxies should experience consider- 
ably more dynamical evolution in a low density 
universe than in the Oq = 1 case. 

A quick estimate whether tidal heating can turn 
a disk galaxy into an elliptical is the increase of the 
velocity dispersion induced by the energy change 
(eq. [18]): = {AE). For a typical tidal 

parameter /tid = 10^ Gyr~^, the inner region at 3 
kpc will be heated to trtid = 54 km s~^. This is not 
enough to completely randomize stellar velocities 
in a large galaxy, but enough to turn a thin disk 
into a thick one (cf. §7.2). In order to reach the 
velocity dispersion of a large elliptical galaxy, 250 
km s~^, the tidal heating parameter would need to 
be 2.2 X 10^ Gyr~^. This is above the maximum 
for our clusters, and therefore the ellipticals are 
unlikely to be created by tidal interactions. 

7.1. Tidal Effects from Close Encounters 

The smoothing filter used in the calculation 
of the tidal field (Appendix B) automatically re- 
moves the effect of close galactic encounters. Their 
additional contribution to the tidal heating rate 
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can be estimated as follows. The number of close 
encounters per galaxy in the simulation as a func- 
tion of impact parameter R = Rene is (cf. eq. 
[131) 

The amplitude of the tidal force due to a galaxy 
of mass Mg is of the order Fais ~ GMg{R)/R^. 
The duration of the encounter is approximately 
R/Venc: SO that the contribution to the heating 
parameter at a distance R is 



/enc(-R) — 




The total contribution of all encounters is the 
integral over all impact parameters weighted by 
the number of encounters, A^enc,!- The perturbing 
halo can be represented by a truncated isother- 
mal sphere with the density profile given by equa- 
tion (Al). The lower limit of integration is set by 
the condition iVonc.i > 1, i-e. that at least one 
encounter takes place; this gives i?niin ~ 3 kpc. 
The upper limit can be taken as the tidal radius, 
-Rmax ^ 30 kpc, but in the end the heating rate de- 
pends only logarithmically on the integration lim- 
its. 

lenc = / Ienc{R) dNenc,l{R) 

»' -Rmin 

V200kms~V V800kms"V 

where ag is the galactic velocity dispersion. This 
value of /enc is comparable to, but smaller than, 
the characteristic parameter Iq. Thus close en- 
counters contribute an important part but do not 
dominate the tidal heating. 

7.2. Self-consistent TV-body simulations 

These analytical estimates will be substantiated 
in Paper II (Gnedin 2003) which presents self- 
consistent A''-body simulations of several galaxies 
selected from each cluster sample. Each galaxy 
is resimulated with 2 x 10^ particles, supplement- 
ing the internal dynamics with the external tidal 
forces of the cluster. Disks of large spiral galaxies 
survive tidal heating but thicken in the vertical di- 
rection by a factor of two, while their dark matter 



halos are truncated at roughly 30 kpc for ag ~ 250 
km s~^. At the end of the simulation they re- 
semble featureless SO galaxies. In contrast, low 
surface brightness galaxies are almost completely 
disrupted by tidal forces. 

Most of the results can be explained qualita- 
tively by the critical tidal density, corresponding 
to the trace of the tidal tensor via Poisson's equa- 
tion: 47rGyOtid = Faa- It is similar to the instanta- 
neous density of the Roche lobe, but reflects the 
highest tidal peaks along the galactic orbit. This 
density ptid determines the truncation radius of 
the halos and the amount of luminous and dark 
mass stripped by tidal heating. 

Galactic iV-body simulations confirm that the 
amount of tidal heating is set by /tid (or ptid)- The 
evolution of galaxies is stronger on average in Clus- 
ter II. However, reflecting the wide distribution 
of tidal parameters in each cluster (Fig. 17), the 
heating on different orbits within the same clus- 
ter varies as much as between the different cluster 
models. 

8. Conclusions 

Tidal interactions may explain the observed dy- 
namical and morphological evolution of galaxies 
in clusters. Using the simulations of hierarchical 
cluster formation, I calculate the tidal forces along 
the galactic trajectories starting at redshift Zg = 5. 
The tidal field shows a lot of temporal variation. 
-1 The peaks of the tidal force do not always corre- 
GjJ^^)with the closest distance to the cluster center 
and are likely to be due to the infalling groups of 
galaxies or individual massive galaxies. The tidal 
peaks are stronger and more impulsive in the low 

models. Also, in those models the cluster forms 
at earlier epochs and the galaxies experience more 
interaction. Thus overall tidal effects are stronger 
in the low f2o clusters. 

Massive galaxies are strongly clustered at the 
end of the simulations. Together with the remain- 
ing substructure, this significantly increases the 
probability of galactic collisions. The number of 
encounters within 10 kpc is about 10 per galaxy in 
the Hubble time. About 0.1% of such encounters 
is expected to result in a merger, defined by the 
analytical merger criteria in §A. This gives a high 
merger probability of the identified halos, 22-29% 
per halo, although most mergers take place before 
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the clusters virialize. 

The amount of dynamical evolution of galaxies 
induced by tidal interactions can be estimated us- 
ing the tidal heating parameter, Itid- The typical 
values of /tid ~ 10"^ Gyr^^ result in an increase 
of the stellar velocity dispersion by ~ 50 km s~^. 
At present this is insufficient to turn a large disk 
galaxy into an elliptical galaxy. However, in the 
future a further secular evolution and transforma- 
tion into the ellipticals seems inevitable. Close 
galactic encounters cannot be resolved in present 
simulations, but they are estimated to contribute 
a fraction of 10% to 50% of the total heating rate. 
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A. Merger Criteria cind Merging Time 

This appendix describes the merger criteria used in §5.3. The physical processes leading to merger are 
the gravitational capture of two galaxies and the subsequent dynamical friction of their halos. The energy 
exchange during the first encounter determines whether the galaxies form a bound pair, while dynamical 
friction determines the timescale for final coalescence. Also, the tidal forces of the cluster can disrupt even 
a bound system at its point of maximum separation. 

For clarity, I model the halos as truncated isothermal spheres with the parameters scaled to the initial 
galaxy mass, Mj. Assume that the velocity dispersion scales as ag = a^{Mi/M^Y/^ , in agreement with the 
Tully-Fisher and Faber-Jackson relations. Here a* = 250 km s~^ and M* = 2 x 10^^ Mq are the fiducial 
circular velocity and mass of an L* galaxy. Furthermore, assume that the velocity dispersion does not 
change in the course of evolution but the density distribution is truncated at a radius Rt = Rt,* 'Jgl'y*- The 
latter assumption is confirmed by the self-consistent A''-body simulations presented in Paper II, and is 
determined from the simulations. The density profile of the halos is 

= ^i^G$Twr '''' ' < ^^'^ 

where Rc is the core radius. The enclosed mass scales almost linearly with Rt . Let galaxy 1 be more massive 
initially than galaxy 2; it will still be more massive after the truncation, Mi > M2- As usual, the two-body 
problem can be written in terms of the reduced particle with the mass /x = M1M2 / {Mi + M2) ■ The dynamics 
is characterized by the radius of encounter, Rene — ^i — r2, and the relative velocity, Venc = vi — V2. 

A.l. Orbital Energy 

Dynamical friction works when the halos overlap. For that the distance of closest approach during the 
first encounter, iienc, must be smaller than the tidal radius of the larger halo, Rt,i. Then the potential 
energy of the binary system is given by 

W = [ $i(r + Re„c)p2(r)dr = -fl^I, (A2) 



where X2 = 72arctan72 ^, 72 = Rc,2/Rt,2 and 



1= - ^5 n / sin6'd6' 



Rc,i ^ |r + Re„cl , 1, + 
1 — zr r arctan 1 — m ■ 



|r-hRenc| Rc,l 2 (r -h Renc)2 + i?c 1 



(A3) 



where x = r/i?t,2- The integral I varies slowly with Rene and is always of order unity. In the limit 

Rc.i,Rt,2 ^ Rt,i and also in the case Rt^i Ri Rt,2, the integral evaluates to X Ri (1 — Xi)(l — X2)- Thus the 
orbital energy of the reduced particle is 

Eorb « Vl, -M2al{l-xi). (A4) 

Note that it does not depend on the radius of encounter as long as the two halos overlap. Also, being enclosed 

within the larger halo, the smaller halo is further truncated according to i?t 2 ~ -Rone <^2/<^i- 

The dissipation of orbital energy into the internal stellar motion can produce a bound orbit even if initially 
Eoih > 0. This is a non-linear interaction and cannot be calculated analytically. Instead, an approximate 
fimctional dependence can be derived and then fitted to the results of A'^body experiments. The velocity 
drag during the encounter follows from the dynamical friction analysis (Binney & Tremaine 1987; also S. 
Tremaine, 1998 private communication): 

dVenc 47rG2M2pi 
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AE ~ MV^ei -F At ^ M it/ ttT^- (A6) 



where the slow- varying Coulomb logarithm, is neglected. The encounter lasts of the order At ^ Renc/Vei 
so that the energy change is 

The last factor assures that equation (A6) remains valid even in the limit where the relative velocity is lower 
than the velocity dispersion of the larger galaxy, Venc ^ ci . 

The normalization is provided by the simulations of Makino & Hut (1997) who considered a hyperbolic 
encounter of two equal-mass galaxies. I derive AE as the initial orbital energy leading to a bound pair of 
galaxies after the first passage. For the encounters with the impact parameter close to the galactic tidal 
radius, the correction factor for equation (A6) is about 2.4. Thus the galaxies become bound if they satisiy 
the capture criterion: 

A. 2. Dynamical Friction Time 

Even if two galaxies satisfy the capture criterion, they may not have enough time to merge until the 
present. To estimate the merger timescale via dynamical friction, assume that the velocity of the reduced 
particle is w « {af + tI)^/^. For circular orbits of radius r, the rate of loss of the specific angular momentum 
is given by the Chandrasekhar formula (eq. [A5]): 

dJ dv M2 a\ 

-dt^'Tt^-M[V-- 

For eccentric orbits, r should be taken as the distance of closest approach. Recent work extends Chan- 
drasekhar's calculation and includes a time-dependent response to the perturbation (Seguin & Dupraz 1994; 
Colpi & Pallavicini 1998; Dominguez-Tenreiro & Gomez-Flechoso 1998). However, the complicated dynamics 
of the interaction cannot be easily parametrized. Instead, I employ the simple functional form given by equa- 
tion (A8) and use numerical simulations again to obtain correct normalization. Taking now dJ/dt « v dr/dt, 
yields the inspiraling time 



tDF — . (A9) 

cri \ 0-2 / 

The normalization follows from the simulations of Dubinski, Mihos & Hernquist (1999) who studied the 
formation of tidal tails by the interacting equal-mass galaxies. From the merging time of their models (their 
Table 3), I deduce top oc r and t^F 3c af""'', roughly in agreement with equation (A9). The resulting 
timescale is 

\ -1 1 /I , 2\ 3/2 

10 kpcy V25O km s"V 

where rj = a2 /(Ji ■ If the time of collision tcoU does not satisiy the dynamical friction criterion tcon+to f <tH, 
the merger is disallowed. 



tDF = 0.6{j^^]{ „i ) -3 Gyr, (AlO) 



A.3. Tidal Disruption 

According to numerical simulations of isolated mergers (Mihos 1999), after a first encounter the galaxies 

separate on a highly elongated orbit before coming back for a final unification. In the cluster, the external 
tidal force at the distance of maximum separation can prevail the gravitational attraction of the two galaxies 
and disrupt the pair. 

Tidal disruption in clusters is important for the following reason. Halos are truncated at the radius Rt^i 
where the tidal force balances the gravitational pull of the galaxy. If after the encounter the smaller galaxy 
moves away more than i?( 1 from the center of the larger galaxy, the external tidal field would be stronger 
than the mutual attraction. 
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There is still a possibility, however, that the pair will merge. The tidal field varies along the galaxy 
trajectory and the truncation is determined by the peak values. If the collision takes place when the external 
tidal field is weak, the galaxies may remain bound. 

Quantitatively, the distance of maximum separation can be determined from the orbital energy, which 
must now be negative after including the dissipated term (eq. [A7]). If the halos detach completely after 
the encounter, 



Eor. = — ^ < 0. (AH) 



The angular momentum term can usually be ignored, so that the maximum distance is 

^ « (1 - xi) (1 + M2/M1) (A12) 

If the value of rm determined from equation (A12) is smaller than Rt,i, equation (All) does not apply and 
the galaxies never separate but undergo a quick merger. 

The external tidal force on the larger galaxy can be written as i^tid = f fmO'i I Rt ii where / = 1 corre- 
sponds to the peak value responsible for the truncation. The galaxy pair remains bound at the distance of 
maximum separation if it satisfies the binding criterion: 



(1-Xi)(l + M2/Mi) 



> 1. (A13) 



B. Calculating the Tidal Field 

This appendix describes the procedure for calculating the tidal field in the simulations, as well as some 
numerical tests. The second derivative of the potential is evaluated using tri-linear interpolation from the 
eight nearest mesh points enclosing the center of mass of the galaxy. In order to remove the effects of small- 
scale noise, I apply the Savitzky-Golay smoothing filters of the 4th order (Press et al. 1992, §14.8, routine 
savgol). These filters replace each data point $i with a linear combination of itself and several neighbors: 

Nr 

J2 <^^n^^i+n, (Bl) 

n=-NL 

where A''^ and are the number of points used "to the left" and "to the right" , respectively. The first and 
second derivatives of the potential are obtained through similar expressions with different coefficients c^^ 
and The double derivatives, F^x, Fyy, F^z, are calculated at the mesh points {i,j,k) as 

, Nr 

F,,{i,j,k) = ^ J2 c(f)$(i + n,i,fc), (B2) 

n=-NL 

where A is the mesh size. The mixed derivatives, F^y, F^z, and Fyz, are calculated using the smoothing in 
both directions: 

-. Nr Nr 

F.y{hhk) = -^ J2 E cWc«$(i + ni,j+n2,fc). (B3) 

ni=—NL Ti2 =— JVi 

The coeflacients are chosen such that the 4th derivative (the smoothing order) of the field is preserved. I take 
Nl = Nn = 4, which corresponds to the effective smoothing of the tidal field on the scale of 250 — 500 
kpc. 
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The accuracy of the above smoothing has been tested against the analytic NFW model (Navarro, Frenk 
& White 1997) with the parameters of the Cluster II simulation: 



p(r) = 



(r/r,)(l + r/r,)2- 



(B4) 



For the virial radius it!vir ~ 1-5 Mpc and the "concentration" relation i?vir = 10 rg, the scale radius is 
Ts = 0.15 Mpc. In the 512^ simulation the grid cell size is A = 62.5 kpc, so that ~ 2.5 A. The 
exact tidal force in the NFW model is readily calculated using equation (5) of Gnedin, Hernquist & Ostriker 
(1999). 

Figure 18 shows that our method reproduces well the mixed derivative Fxy down to 3 cell sizes from the 
center, and somewhat worse the double derivative F,,;,,.. In the case; of the mixed derivative, the filter may 
seem to be less accurate than the straightforward two-point finite difference. But in real simulations with 
numerical noise the smoothing filter becomes more effective and removes small-scale artifacts. 

B.l. Test: The Pancake Simulation 

As another check of accuracy of the calculation of the tidal field, I have run a test simulation of a one- 
dimensional Zeldovich's pancake. The initial conditions are set up as described in Bryan et al. (1995). The 
wavelength of the sine wave perturbation equals the size of the simulation box; the maximum displacement 
is in the center. The simulation is run until a caustic of infinite density forms. From Poisson's equation 
it follows that the only non-zero component of the tidal tensor is F^x = 4:-jTGAp, where Ap is the density 
perturbation. 

I have run two simulations, with the high (512^ grid) and low (128^ grid) resolution. The density profile 
in both simulations is reproduced accurately. The calculated tidal force matches the exact solution in the 
beginning while the perturbation is small, but deviates at later times. The smoothing filter cannot resolve 
a very narrow caustic. As expected, the higher-resolution simulation departs from the exact solution later 
than the lower-resolution one, at 70% vs 20% of the total simulation time. A similar effect of the resolution 
loss in high density regions in real simulations is emphasized again in the next section. 

B.2. Subtracting Galactic Self-Contribution 

Large halos may contribute significantly to the amplitude of the local tidal field. If the halo extends over 
more than a couple of cells, it will contaminate the calculation of the external tidal force. It is therefore 
necessary to correct for the galactic self-contribution. 

For this purpose, I have implemented a potential solver for an isolated galaxy. On a rectangular grid of 16^ 
cells, the density of the galactic particles is assigned using the CIC method. Then, the potential is obtained 
by the FFT method with isolated boundary conditions. Applying the same procedure for calculating the 
tidal tensor as above, I evaluate and subtract the halo self-contribution. 

The FFT solver for Poisson's equation with isolated boundary conditions is implemented following Hock- 
ney & Eastwood (1981). The density grid is doubled in all three directions and padded with zeros in the 
added space. Green's function is translated correctly to the new sections in real space and then Fourier 
transformed and stored for the potential calculation. The trick with extending the grid is necessary in order 
to avoid the multiple images which appear in the usual FFT with periodic boundary conditions. The result- 
ing Green's function Gisoi{k) differs from the usual form, G{k) = k~^, by cutting off the high frequency tail. 
Afterwards the algorithm proceeds as usual: Fourier transforming the density field, multiplying by Giso\{k), 
and transforming back to real space. 

The isolated FFT procedure has been tested on the same 16"^ grid, using 2^^ particles distributed in the 
central cell according to the Hernquist (1990) model: 




Ma 1 



27r r {r + a) 



,3' 



(B5) 
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with a tiny core radius, a = 0.1 A. Essentially, this is a point mass represented by particles. The algorithm 
reproduces the exact potential with very high accuracy, all the way to the edge of the simulation box. (For 
comparison, the periodic FFT produces large errors on most of the grid.) The force calculation is limited 
in the center at 2 A, as in any grid-based method, but behaves perfectly well at larger radii. The force 
calculated with the periodic FFT shows much more noise. 

Another test has been applied to make sure the procedure works in a real simulation. The 16^ density 
field was calculated using the particles from a large galaxy, which were all located within a single grid cell 
and resembled a smoothed point source. Although with increased scatter, the new method is still in a good 
agreement with the exact potential and force. 

Figure 19 shows the self-contribution of a large galaxy in the Cluster III simulation. For the ease of 
comparison, I plot the trace of the tidal tensor (Fxx + Fyy + F^z) which is generally proportional to the 
density of galactic particles, in accord with Poisson's equation. The smoothing filter misses the extreme 
variations of the density, but otherwise the amplitude of the tidal force is calculated correctly. 

Also, Figure 19 shows the scaling (1 + z)'^, which corresponds to a constant comoving density. The halo 
density follows this scaling quite closely indicating that the galaxy expands in real space. This numerical 
artifact is due to the loss of force resolution within a grid cell. The particles cannot collapse below 1 A, and 
therefore in comoving coordinates the galaxy continues to occupy the same volume. 
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Fig. 1. — The surface density map of the three 
simulated clusters at ^; = 0. The projection is 
over the X direction. Saturated white regions are 
denser than 200 times the mean density; regions 
below the mean density are black. The size of the 
windows is 8/i~^ Mpc. Thick bars in the lower 
left corner of each panel indicate 1 Mpc. The 
lower right panel shows the enlarged view of the 
inner region of the VIq — 1 cluster. 



This 2-column preprint was prepared with the AAS lAlJ^X 
macros v5.0. 
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Fig. 2. — The density profile of the three clusters. 
Dashes and open squares are for the Oq = 1 model, 
dots and filled triangles are for the = 0.4 
model, solid line and stars are for the — 0.4, 
f^A = 0.6 model. The innermost point corresponds 
to the resolution element, A. 
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Fig. 4. — Resolution effects on the structure of the 
Clo = 0.4, Qa = 0.6 cluster. Solid Hue and stars 
are for the main 512"^ grid cells simulation, dashes 
and open squares arc for the 256'^ simulation, and 
dots and filled triangles are for the 128"^ simula- 
tion. Arrows point at the respective cell sizes. A, 
in each simulation. The nominal PM force resolu- 
tion is 2A. 



Fig. 3. — The velocity dispersion of dark matter 
in the three clusters. Line notation as in Figure 
2. Note a rise of the velocity dispersion at 2 
Mpc from the center of the = 1 cluster due to 
the new infalling group of galaxies. 
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Fig. 5. — Three-dimensional distribution of parti- 
cles in the identified halos at Zg = 5 in the fJo = 1 
cluster. The 100 most massive halos, whose tra- 
jectories are followed in the simulation, are black; 
others are gray. Boundary size is 16 Mpc co- 
moving. 




Fig. 6. — Particles in the identified halos in the 
flo = 0.4 cluster. 



Fig. 7— Particles in the identified halos in the 
Qo = 0.4, riA = 0.6 cluster. 
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Fig. 8. — The mass function of identified halos 
in the three clusters (histograms). Superimposed 
axe the predictions of the Press-Schechter theory 
for the corresponding power spectra (lines). The 
lower right panel shows, for comparison, the mass 
function from the lower resolution, 256^ simula- 
tion of Cluster III (shaded histogram). 
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Fig. 9. — Three projections of galaxy trajectories 
in the Oq = 1 cluster (comoving coordinates). One 
galaxy is very massive (solid line), ranked ^4 by 
mass in the list of all identified halos, the other is 
ranked #93 (dotted line). The lower right panel 
shows the distance from the final cluster center. 
Large dots mark the initial positions where the 
halos were identified. 
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Fig. 11. — Three projections of galaxy trajectories 

in the f^o = 0.4, f^A = 0.6 cluster. Notation as in 
Figure 9. Solid line is for the galaxy ranked #4, 
dotted line is for the galaxy ranked #87. 
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Fig. 10. — Three projections of galaxy trajectories 
in the VLq = 0.4 cluster. Notation as in Figure 9. 
Solid line is for the galaxy ranked #4, dotted line 
is for the galaxy ranked #93. 
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Fig. 12. — Distribution of close encounters in time: 

Number of encounters with i?,cnc < 20 kpc per 
galaxy per Gyr. Filled histograms show the num- 
ber of possible mergers. 
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Fig. 13. — Distribution of relative velocities: 
Number of encounters with i?enc < 20 kpc per 
galaxy per 100 km s~^. 



Fig. 15. — Tidal field around the two galaxies in 
the Clo = 0.4 cluster (cf. Fig. 10). Line notation 
as in Figure 14. 
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Fig. 14. — Tidal field around the large (solid line) 

and dwarf (dots) galaxies in the f^o = 1 cluster 
(cf. Fig. 9). Upper panel shows the distance from 
the final cluster center. 
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Fig. 16. — Tidal field around the two galaxies in 
the Qo = 0.4, Oa = 0.6 cluster (cf. Fig. 11). Line 
notation as in Figure 14. 
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Fig. 17. — The distribution of the tidal heating 
parameter, /tid (eq. [19]), for the identified galax- 
ies in Cluster I (filled histogram) , Cluster II (solid 
histogram), and Cluster III (hatched histogram). 




5 10 15 



r/A 

Fig. 18. — The tidal tensor components in the 
NFW model (arbitrary units) versus radius in 
units of the cell size. Dashed line is obtained 
with the direct two-point finite difference, solid 
line with our smoothing filter. Dots show the an- 
alytic solution. Note that F^x changes sign at 
r « 5 A; the inner curve shows the absolute value. 
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Fig. 19. Galaxy self-contribution to the trace 
of the tidal tensor (solid line) compared to the 
density of galactic particles {AnGp; dots) of halo 
#4 in the Cluster III simulation. Dashes show 
the (1 + z)^ scaling for a fixed comoving density 
(arbitrary normalization). 
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Table 1 
cosmological models 



Cluster rio h 



I 1 0.65 0.5 

II 0.4 0.65 1.0 

III 0.4 0.6 0.65 1.0 



Table 2 

Parameters of the simulated clusters 











b/a 


c/a 


Cluster 


(/i-i Mpc) 




(km s~^) 






I 


1.03 


3.9 X lO^"* 


656 


0.70 


0.58 


II 


1.47 


4.5 X lO^"* 


664 


0.85 


0.71 


III 


1.42 


4.1 X 10^* 


649 


0.85 


0.73 



Table 3 

Resolution test on the parameters of Cluster III 







Mvir 




b/a 


c/a 




{h-^ Mpc) 


(Mo) 


(km s-i) 






512^ 


1.42 


4.1 X 10" 


649 


0.85 


0.73 


256-'^ 


1.42 


4.1 X lO^'' 


632 


0.84 


0.70 


128^ 


1.43 


4.1 X 10" 


546 


0.76 


0.61 



Table 4 

Parameters of calaxy distribution 



Cluster 




ro {h ^ kpc) 


7 


C 


All galaxies 


I 


34 


11 


2.4 


68 


II 


69 


70 


4.2 


510 


III 


65 


22 


3.4 


1600 


Galaxies with r > A = Q2.5 h ^ 


kpc 


I 


33 


66 


2.9 


15 


II 


61 


110 


4.5 


230 


III 


4:-! 


.19 




100 



28 



Table 5 

Number of mergers and close encounters 



Rene < 5 kpc 10 kpc 15 kpc 20 kpc Total 



Cluster I 



Nenc 


32 


97 


139 


211 


498 







5 


8 


8 


10 


Nenc,l 


1.9 


5.7 


8.2 


12 


29 







0.15 


0.24 


0.24 


0.29 


^enc 


32 


130 


241 


381 


890 


A^enc (test) 


3 


23 


43 


80 


213 


-/Vmer (test) 














1 


Cluster II 


Nenc 


435 


866 


1128 


1280 


2576 


mer 





9 


18 


20 


20 


-^enc,l 


13 


26 


33 


38 


76 


mer,l 





0.13 


0.26 


0.29 


0.29 


^cnc 


435 


1404 


2596 


3590 


6631 


A^enc (test) 


19 


62 


143 


273 


902 


-/Vmer (test) 

















Cluster III 


Nenc 


95 


314 


575 


782 


1877 




4 


7 


10 


11 


14 




3.0 


9.8 


18 


24 


59 




0.06 


0.11 


0.15 


0.17 


0.22 


^ *enc 


136 


525 


962 


1360 


3113 


Nenc (test) 


13 


59 


129 


197 


614 


A^mer (test) 














1 



Table 6 

Relative velocities in close encounters 







Kied 




Cluster 


(km s-i) 


(km s-1) 


(km s-i) 


I 


550 


686 


936 


II 


350 


872 


1196 


III 


350 


817 


1097 



29 



Table 7 
Galaxy Kinematics 



Region < r Ng{r) cjg^ios 
{h-^ Mpc) (km s-i) 



Cluster I 



Inner 


0.21 11 


373 


Middle 


0.65 11 


324 


Rvir 


1.03 12 


538 


Cluster II 


Inner 


0.11 23 


117 


Middle 


0.34 23 


433 


Ry'lT 


1.47 23 


469 


Cluster III 


Inner 


0.063 22 


131 


Middle 


0.5 22 


324 


Ry'lT 


1.42 21 


517 



30 



